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Abstract 

We report results of statistical and dynamic analysis of the serrated stress-time curves obtained from compressive 
constant strain-rate tests on two metallic glass samples with different ductility levels in an effort to extract hidden 
information in the seemingly irregular serrations. Two distinct types of dynamics are detected in these two alloy 
samples. The stress-strain curve corresponding to the less ductile Zr^ Cu 15M10A/10 alloy is shown to exhibit finite 
correlation dimension and a positive Lyapunov exponent, suggesting that the underlying dynamics is chaotic. In 
contrast, for the more ductile CM47.5Zr47.5AZ5 alloy, the distributions of stress drop magnitudes and their time durations 
obey a power law scaling reminiscent of a self-organized critical state. The exponents also satisfy the scaling relation 
compatible with self-organized criticality. Possible physical mechanisms contributing to the two distinct dynamic 
regimes are discussed by drawing on the analogy with the serrated yielding of crystalline samples. The analysis, 
together with some physical reasoning, suggests that plasticity in the less ductile sample can be attributed to stick-slip 
of single shear band, while that of the more ductile sample could be attributed to the simultaneous nucleation of large 
number of shear bands and their mutual interactions. 
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1. Introduction 

The discontinuous yielding of crystalline metallic al- 
loys during plastic deformation is well known to exhibit 
temporal fluctuations in the stress-strain curves. These 
fluctuations are accompanied by spatial localization of 
dislocations into bands U]. This phenomenon, known 
as the Portevin-Le Chatelier (PLC) effect [2], has been 
studied extensively in crystalline materials, and contin- 
ues to attract considerable attention lHH El S H H 
[loll . The physical origin of the intermittent flow is at- 
tributed to the collective behavior of dislocations |ll|3|. 
In contrast, the less-well known, but similar serrated 
flow along with localization of deformation into shear 
bands has been observed in a large number of stud- 
ies on bulk metallic glasses (BMGs) despite the fact 
that BMGs are characterized by the absence of long- 
range order. Although dislocation-mediated deforma- 
tion mechanism is absent in BMGs, serrated flow phe- 
nomenon in the plastic regime has been widely observed 
in nanoindentation experiments 111 l[ll2[ll3[ll4ll5ll and 



compression tests H [nl El H I2III . The plas- 
tic deformation of BMGs at room temperature is char- 
acterized by nanoscale shear bands that contribute to 
strain softening arising from dilation [22], or an in- 
crease in temperature, or both [17l 23 1. These serra- 
tions have been attributed to shear-banding operation. 
As the strain rate increases, the amplitude of serrations 
decreases with fewer shear bands simultaneously con- 
tributing to serrations at high strain rates compared to 
their number at low strain rates N19II . In the case of 
nanoindentations, the serrations have been reported to 
disappear at high rates of indentation. This has been 
attributed to the inability of a single shear band to ac- 
commodate the entire shear at high rates 01 IB. How- 
ever, Jiang and Atzmon IU3I1 and Greer et al. 11411 have 
considered the disappearance of serrations at high strain 
rates as arising from limitations of data acquisition or 
due to proliferation of shear bands. Subsequent studies 
demonstrated that the deformation remains inhomoge- 
neous even at high indentation rates 11511 . which is also 
consistent with results on compression tests [19]. 
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The inhomogeneous deformation of metallic glasses 
has also been studied theoretically. The conventional 
picture of the deformation of amorphous materials is 
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based on the free volume model |22, 24ll o r the shear 
transformation zones (STZ) model |25, 261. which as- 
cribe the macroscopic deformation to the presence of 
localized deformation centers. This is also supported 



by numerical simulations [26, 27, 28, 29, 30, 31]. Such 



shear transformations should lead to long-range inter- 
action. Models that include such interactions do exhibit 
localized deformation patterns that resemble the shear 
bands I29il32ll33ll . Another direction of investigation of 
intermittent plastic flow recognizes the importance of 
nonlinearities that naturally lead to loss of stability of 
the "elastic branch" or equivalently to loss of analytic - 
ity of stress [34]. 

Here we adopt a completely new approach to analyze 
the serrated flow of BMGs in an effort to obtain insight 
in the deformation dynamics of these materials. In gen- 
eral, the serrated stress-strain curves of BMGs appear 
quite complex and irregular. The question we address 
in this paper is whether the seemingly irregular nature 
of stress-strain curves is purely random or there is any 
hidden order that could extracted from them and what 
physical conclusions could be drawn from such a study. 
Indeed, such efforts have been carried out on serrated 
flow of the crystalline alloys, namely the PLC effect 



nonlinear dynamics, called nonlinear time series analy- 
sis, is designed to recover hidden information from the 
seemingly irregular stress-strain series alone. The well- 
developed subject of nonlinear dynamics suggests that 
any irregular time series could well arise due to non- 
linear deterministic interaction of at least three degrees 
of freedom, referred to aschaos, instead of being purely 
random in origin. Indeed, such irregular time series are 
found in a large number of physical situations including 
physical, chemical and biological systems which have 
been demonstrated to be chaotic 113611 . 

That such an approach could be useful is also sug- 
gested by the fact that the intermittent flow of BMGs has 
been recognized to be a stick-slip phenomenon B371I38I1 . 
Indeed, a large number of systems ranging from atomic 
scale (e.g. atomic force microscopy) to the geological 
scale (earthquakes) exhibit stick-slip dynamics. At lab- 
oratory scale, stick-slip phenomenon has been observed 
in a large number of rate-controlled experiments such 
as the PLC effect I3l 14. 13911 . peeling of an adhesive tape 
||40L 41], micellar vesicles subjected to shear 042ll . etc. 
The phenomenon is characterized by a sticking phase 
that usually lasts much longer than the slipping phase, 
a feature also observed in the serrations of BMG defor- 
mation. Stick-slip dynamics is well established to re- 
sult from a competition between the internal relaxation 
time scales and the time scale associated with drive rate 



1 4311 . Indeed, this is consistent with the fact that sev- 
eral time scales corresponding to free volume, temper- 
ature, etc., are well known to influence BMG deforma- 
tion II22L 12511 . While the above facts suggest that such 
an analysis is worthwhile, further support also comes 
from the recently established similarity between the ser- 
rated flow in BMG deformation and the PLC effect for 
crystalline alloys 04411 . The latter is a well established 
example of stick-slip dynamics arising from the collec- 
tive pinning and unpinning of dislocations from solute 
atmosphere QjJ, |3j |4J] . 

In view of the similarities between the serrated flow 
in BMGs and the PLC effect, we recall some rele- 
vant methods adopted in analyzing the serrated flow 
of the PLC effect. Time-series analysis of the stress- 
strain curves from the PLC effect of single crystals and 
polycrystals have been carried out. Such an analy- 
sis was motivated by the prediction by the Ananthakr- 
ishna model (1 0, IE 51 III H that the stress-strain 



curves can be chaotic for a range of strain rates£] The 
results of the analysis show the existence of two dis- 
tinct types of dynamic regimes, namely chaotic and 
self-organized criticality [8, 9, 1( 



3511 . in two different 



35J]. The methodology rooted in the area of tion O, |4j, |fl |8j, |9l [1 



regimes of strain rate, consistent with the model predic- 
45, H 47]. So far such a system- 



atic dynamic analysis of the serrated flow of the BMGs 
has not been performed. It is the purpose of this paper 
to carry out dynamic time-series analysis of the stress- 
strain curves of two different types of BMGs. We closely 
follow the methods used in Refs. fl H, OH B . Si- 
multaneously, we have carried out a detailed statistical 
analysis of stress drop magnitudes and durations iKstl . 
Recently Wang et ai, showed that the statistics of the 
elastic energy density of several BMGs J48ll could be 
fitted to a power-law distribution modulated by a Gaus- 
sian function, suggesting that the BMG is in a state of 
self-organized criticality (SOC) 1]. However, 

Wang et ai, did not carry out a full analysis involv- 
ing the distributions of the durations of the events, as 
well as the consistency checks essential to assert SOC 
dynamics. Our detailed analysis shows that while the 
less ductile Zr^Cu^NixoAlxo BMG is chaotic, the rela- 
tively more ductile sample Cu^^Zr^ 5A/5 surprisingly 
exhibits a power-law state of stress drops, suggesting 
self-organized critical state. Physical mechanisms con- 
tributing to these two distinct dynamic regimes are dis- 



Indeed, the Ananthakrishna model uses a dynamic approach to 
study the PLC effect 1451 |46L 14711 and considers the phenomenon 
as arising from nonlinear interaction of a few collective modes of 
dislocations 14611 - The model recovers most generic features of 
the PLC effect apart from predicting the serrations to be chaotic 
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cussed. Our analysis also suggests that the high ductility 
of Cun ,sZr47 5A/5 sample emerges from the production 
of several shear bands and their interactions rather than 
partial plastic relaxation of the shear activity as in the 
PLC effect fl, Si 



2. Experimental procedure 

Zr^Cu\sNi\oAlu) and CM47.5Zr47.5AZ5 alloys with 
nominal compositions were prepared by arc melting a 
mixture of pure metal elements in an argon atmosphere, 
followed by suction casting into copper molds to form 
rod-like BMG samples with a diameter of 2mm. The 
glassy phase of the as-cast BMG specimens were exam- 
ined by transmission electron microscope (TEM) and 
X-ray diffraction (XRD). The specimen surface was ex- 
amined by scanning electron microscope (SEM; Philips 
XL 30). The compressive test specimens were cut from 
the BMG rods with a nominal 2:1 aspect ratio. The 
samples were carefully ground to Lum accuracy. Com- 
pressive tests were performed in a Instron electrome- 
chanical 3384 test system at a constant strain rate of 
1 x 10~ 4 s~ l . Care was exercised to ensure that the two 
surfaces were nearly parallel and orthogonal to the load- 
ing axis. Load, displacement and time were measured 
and digitally stored at a frequency of 50 Hz (50 datum 
point per second) for the Zr(,sCu\sNin)Al\Q BMG sam- 
ple and 2QHz for the CM47.5Zr47.5AZ5 BMG sample. 
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Figure 1: Stress-strain curves for the two BMGs of composition 
Zr^Cu^NiioAlio and Cu^^Zr^^Al;, respectively. 

Figure[T]shows stress-strain curves for the two BMG 
samples. The lower curve corresponds to the sample of 
composition ZrgsCKisMioA/io and the upper curve cor- 
responds to the sample of composition CM47.5Zr47.5AZ5. 
It is clear from the figure that the two alloys dis- 
play significantly different plastic strains (about 8% for 
ZresCMisMioAZio and 11% for CM47.5Zr47.5AZ5) before 
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Figure 2: (a,b) A small portion of the stress-strain curves for the 
Zr^sCuisNiioAlio sample and CH47 sZr^ 5A/5 sample, respectively. 



the catastrophic failure. Thus, for convenience on no- 
tation, we label the less ductile Zr^CuisNimAlio BMG 
sample by S b and we shall generally refer to the more 
ductile CM47.5Zr47.5AZ5 BMG sample S d- Clearly, both 
samples exhibit intermittent flow after yielding at about 
2% elastic strain. The serrations are characterized by 
the repetitive cycles of abrupt stress drops followed by 
elastic reloading (Fig. [2^ and b). As can be seen, the 
serrations for the Sb (Zr^sCuisNiioAlio) sample have 
large stress drops with very few small ones, while the 
serrations of the more ductile sample So exhibits nu- 
merous small stress drops with fewer large ones. Plots 
of stress-strain curves of the two BMGs for the same 
strain range are shown in Figure |2] We note that sam- 
pling rate is such that for both samples, a few points are 
recorded on the region of stress drops. 

Figure [3] shows the shear band patterns on the surface 
of the deformed samples for the two BMGs. As shown 
in Fig. [3^, a dominant shear band was observed dur- 
ing the deformation of the less ductile Zr^Cux^NixoAlxQ 
BMG sample, even though a fair amount of ductility is 
exhibited by the sample. The reasonably large plasticity 
when a single dominant shear band operates has been 
recently attributed to stick-slip process 
contrast, multiple shear bands were seen on the sur- 
face of the deformed sample for the more ductile BMG 
CM47.5Zr47.5AZ5 as shown in Fig. [3}}. A magnified image 
that displays numerous shear bands when the sample is 
deformed significantly is shown in Fig. [3};. The simul- 
taneous nucleation of large number of shear bands is 
responsible for the good plasticity of the ductile sample 
Cw47.5Zr47.5AZ5 II53I1 . (See section 5 for further discus- 
sion.) 

In general, work or strain hardening is not observed 
during the deformation of the BMGs as shear bands act 
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Figure 3: (a and b) Shear bands for the less ductile Zr^sCuisNi^Alio 
BMG sample and the more ductile Cut,-] jZr^ 5A/5 sample, (c) Mag- 
nified view of the more ductile metallic glass where intersection, ar- 
rests and branching of the shear bands are shown. 



as work-softening sources attributed to shear dilation or 
temperature rise. In the present case, however, a clear 
upward drift in the stress-strain curve of the ductile sam- 
ple S d can be observed. This "work hardening" behav- 
ior is expected to originate from the extensive multi- 
plication of shear bands and consequent interaction be- 
tween them ll53ll . Much like hardening in the PLC ef- 
fect, the gradual increase in the stress level is also ac- 
companied by increase in the magnitude of stress drops 
with increasing strain. 

3. Methods of analysis 

Given a time series, two kinds of analyses are possi- 
ble. The first one is the statistical analysis of the time 
series. While the method is simple, it can yield substan- 
tial information on the nature of the underlying process. 
The second method is the more sophisticated method of 
nonlinear time-series analysis. This approach is particu- 
larly useful if one suspects that the underlying process is 
of deterministic origin. As stick-slip process is a deter- 
ministic nonlinear phenomenon, the method of nonlin- 
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Figure 4: (a and b) Cleaned stress-time series for the less ductile 
Zr^jCu[^Ni[oAl[o BMG and the more ductile CM47 ^Zr^-j 5A/5 BMG 
samples, respectively, for the region indicated in the original stress- 
strain plot. 



ear time-series analysis offers a platform for analyzing 
the phenomenon and also to understand the complexity 
of the stick-slip process. It will also help to establish if 
the irregular time series is chaotic. Both approaches as- 
sume infinitely long time series. We shall comment on 
the applicability of the methods to real systems. 

3.1. Statistical analysis: Power law distributions 

The simplest statistical quantity to compute is the 
statistics of events. The definition of an event depends 
on the physical situation, which in our case may be 
taken to be the magnitude of the stress drop. Since the 
stress drop of a certain magnitude occurs in a certain 
time interval, the duration of the stress drop is also a 
characteristic feature of the event. In our analysis, we 
have taken the difference between the maximum and 
next minimum to be a measure of the magnitude of the 
stress drop denoted by Act. Similarly, we define the time 
interval between a stress minimum to the previous max- 
imum to be the duration of the stress drop denoted by 
Af. It is clear that the corresponding distributions given 
by D(Acr) and D(At) can be easily computed. As is the 
case for the PLC effect HSOU], the distribution of the 
event sizes D(Acr) often follows a power-law scaling of 
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the form 



D(Ao-) ~ Act"", 



(1) 



where a is an exponent. Similarly, D(At) also follows a 
power law scaling of the form 



D(Af) ~ AT 1 *, 



(2) 



with j3 being the corresponding exponent. However, 
since Act refers to the stress drop occurring in a time 
interval Af, these two variables are not independent. In- 
deed, the basic statistical quantity describing the pro- 
cess is the joint probability of an event of magnitude 
Acr occurring in a time interval Af defined by P(A<x, Af) 
05411 . Thus, both P( Acr) and P{At) are the marginal dis- 
tributions obtained by integrating out the complimen- 
tary variable. (Note our D's are unnormalized.) Thus, 
one can introduce another scaling exponent: 



< Af >~ Acr x 



(3) 



where < . > is the conditional average. Then, the three 
exponents are connected through a scaling relation: 



a = x(0- 1) + 1. 



(4) 



3.2. Time series analysis 

Quite often, only a scalar signal is measurable, which 
in our case is the stress. Such signals can be very ir- 
regular. This again is the case in our deformation ex- 
periments. Traditionally, such irregularity is generally 
attributed to statistical noise such as the Brownian noise 
which can only be characterized by their distribution (or 
their moments). However, the well-developed methods 
in dynamic systems suggest that such signals could well 
arise from nonlinear interaction of a few degrees free- 
dom (at least three) called deterministic chaos or simply 
chaos. Two characteristic features of chaos are the ex- 
istence of a strange attractor with self-similar properties 
quantified by a fractal dimension (or equivalently the 
correlation dimension), and sensitivity to initial condi- 
tions quantified by the existence of a positive Lyapunov 
exponent. Given the equations of motion, these quan- 
tities can be directly calculated. In contrast to the sta- 
tistical analysis, time-series analysis is a basic tool for 
establishing if a time series is of dynamical origin or 
not. This method is useful when a scalar time series 
is suspected to be a projection from a higher dimen- 
sional dynamics. Such time series are traditionally ana- 
lyzed using embedding methods that attempt to recover 
the underlying dynamics. The basic idea is to unfold 
the dynamics through a phase space reconstruction of 
the attractor by embedding the time series in a higher - 
dimensional space using a suitable time delay 0551 Dal . 



Consider a scalar time series measured in units of sam- 
pling time St defined by [x(k), k — 1,2, 3, • • -, N]. Then, 
one can construct li-dimensional vectors defined by: 

£k - [x(k), x(k + t), ■ • •, x(k + (d - l)r)], 
k = l,---,[N-(d-l)T]. (5) 

The delay time t suitable for the purpose is either ob- 
tained from the autocorrelation function or from mutual 
information 13611 . Once the reconstructed attractor is ob- 
tained, the existence of converged values of finite corre- 
lation dimension and a positive exponent is taken to be 
a signature of the underlying chaotic dynamics. 

3.3. Correlation dimension 

The self-similar nature of the reconstructed attracted 
is quantified by using the popular algorithm due to 
Grassberger and Procaccia 05611 that calculates the cor- 
relation integral. This is defined as the fraction of pairs 
of points |j and $j whose distance is less than r given 
by: 



p Ui 



(6) 



where ©(■•■) is the step function and N p the number 
of vector pairs summed. A window is imposed to ex- 
clude temporally correlated points B36l 15711 . The method 
provides equivalence between the reconstructed attrac- 
tor and the original attractor. It has been shown that a 
proper equivalence is possible if the time series is noise 
free and long 05 811 . For a self similar attractor C(r) ~ f 
in the limit of small r, where v is the correlation dimen- 
sion l56Tl . Then, as d is increased, one expects to find 
a convergence of the slopes d lnC(r)/dln r to a finite 
value in the limit of small r. However, in practice, the 
scaling regime is found at intermediate length scales due 
to limited length of the time series and the presence of 
noise. 

3.4. Lyapunov spectrum 

The existence of a positive Lyapunov exponent is 
considered as an unambiguous quantifier of chaotic dy- 
namics. However, the presence of superposed noise 
component, which in the present case is high, poses 
problems. In principle, the noise component can be 
cured and then the Lyapunov exponent calculated [36]. 
Here, we use an algorithm that does not require pre- 
processing of the data; it is designed to average out the 
influence of superposed noise. The algorithm, which 
is an extension of the Eckmann's algorithm, has been 
shown to work well for reasonably high levels of noise 
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in model systems as well as for short time series (for de- 
tails, see Ref. lfloll59ll ). The method has also been used 
to analyze experimental time series. 

The conventional Eckmann's algorithm [60] relies on 
the construction of a sequence of tangent matrices T,, 
which maps the difference vector - gj to - 
evolved for k units of time, and successively reorthogo- 
nalize T, using standard QR decomposition (Q is an or- 
thogonal matrix and R is an upper triangle matrix with 
positive diagonal elements). Then, the Lyapunov expo- 
nents are given by: 



1 



.i/ = ^EfrC*,) fl ;Z = i.2.-»* 

7=0 



(7) 



Here p is the number of available matrices and k is the 
propagation time in units of the sampling time 5t. In 
the algorithm, the number of neighbors used is small, 
typically min[2c/, d + 4], contained in a sphere of size 
e s . A simple modification of this is to use those neigh- 
bors falling in shell e s defined by spheres of inner and 
outer radii e, and eo, respectively. Then, ignoring points 
falling within the sphere, e, acts as a noise filter. How- 
ever, these few neighbors will not be adequate to aver- 
age out the noise component superposed on the signal. 
Thus, the modification effected is to allow more neigh- 
bors so that the noise statistics is sampled properly. (For 
details, see Refs. lfloll59ll .) As the sum of the exponents 
should be negative for a dissipative system, we impose 
this as a constraint. In addition, we also demand the 
existence of a stable positive exponent and a zero expo- 
nent (a necessary requirement for continuous time sys- 
tems such as the stress-time series) over a finite range of 
shell sizes e s . (Note that the constancy of the two expo- 
nents in the range of e s is equivalent to constancy with 
respect to the number of neighbors contained in e s as it 
is varied over a certain range.) The results presented be- 
low are for a specific value of the shell size that is in the 
stable range of shell sizes for which positive and zero 
Lyapunov exponents are constant. 

Finally, the Lyapunov dimension can be calculated 
from the Lyapunov spectrum by using the Kaplan- Yorke 
conjecture, defined by the relation: 



2 * 

Dry =J+ - — r, and ) A t > 0. 



(8) 



4. Dynamics of serrated flow in bulk metallic glasses 

From a dynamic point of view, the irregular form 
of the stress-strain curves from the two BMG samples 
could be of deterministic nonlinear origin. This point 
of view is also supported by the similarity of the stress- 
strain curves of the BMGs to those based on the PLC 
effect, coupled with the fact that dynamic analysis of the 
PLC stress-strain curves have been shown to be chaotic 
IH OH B . Further support comes from similar anal- 
ysis of other intermittent flows 1 36 , 40[ 41 , 421. 



4.1. Chaotic dynamics 

A few comments are desirable on the limitations in 
applying the above theoretical approaches developed for 
ideal systems to nonideal experimental stress-strain se- 
ries. Nonlinear time-series analysis assumes that the 
signal is of infinite extent and is stationary. In real- 
ity, however, experimental data are always of limited 
length. More importantly, very often, the stationarity 
assumption is not satisfied, as is the case in the current 
problem. The first restriction is less serious and, indeed, 
algorithms have been developed to deal with short time 
series Hoi l36l 59l 621. Such methods have been ap- 
plied to a number of physical situations ranging from 
physics to biology ID6I1 . The second restriction is more 
serious. Non-stationarity is common in many real sys- 
tems of practical relevance such as stock prices, inter- 
mittent natural phenomenon such as earthquakes, and 
signals from living matter. The specific type of non- 
stationarity relevant for the current study arises from the 
slow variation of internal state variables compared to the 
rapidly fluctuating signal that is being analyzed @. This 
is reflected in the systematic change in the mean stress 
level seen, particularly, in the ductile sample. Neverthe- 
less, if some insight about the underlying physical pro- 
cess is available, which is hardening in our case, this 
contribution can be subtracted to "convert" the non- 
stationary time series into a "stationary one". If the re- 
sulting time series is stationary, various statistical and 
dynamic quantities, such as the mean, variance, more 
importantly the correlation functions and correlation in- 
tegrals should be the same (within some reasonable er- 
rors) for reasonably long subintervals of the whole time 
series. Furthermore, in real systems, most experimental 
signals contain superposed noise, which in this case is 



It has been conjectured that Dky ^ v. For consistency, 
we have also calculated Dky JUl. 



2 A good discussion on tests for stationarity in the context of non- 
linear time-series analysis can be found in Ref. 1 36] . The minimum 
requirement is that the total length of the time series should be much 
larger than some basic time scale, which in this case is the mean pe- 
riod of the stress fluctuations. 
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high. There are several methods designed to cure the 
noise component including singular value decomposi- 



tion [36, 63]. Usually, the cured data sets are then sub- 



jected to further analysis. 

Apart from noise, there are three types of systematic 
variations in the time series. The first one corresponds to 
near-periodic fluctuations in the stress that can be traced 
to the vibrational frequency of the testing machine. The 
second is the accuracy in stress measurement, which is 
less than O.SMPa. The first error is easily eliminated. 
However, the accuracy of measurement is a limitation 
that is reflected in the extent of scaling regime in the 
correlation integral. The third is intrinsic to the material, 
namely the gradual increase (for the ductile sample) in 
the mean stress level or the decrease in the mean stress 
level towards the end of the time series. The gradual 
increase in stress level is known to arise due to harden- 
ing and is considered as not relevant to the mechanisms 
controlling the stress drops. Therefore, this component 
is eliminated by subtracting a running average from the 
original time series. As for the less ductile BMG sample 
S b, the changes in the mean stress level, though small, 
are also dealt with in a similar manner. The stress drop 
magnitudes are also normalized by using a running av- 



erage of the stress drops QSj, |9j [10l |35j . Further, the sam- 
ple eventually breaks, and thus, not only the time se- 
ries is limited, but the eventual failure of the sample can 
induce artefacts such as a decrease in the mean stress 
level. This also implies that we need to discard a certain 
portion of the time series (typically last 10-20 peaks) 
towards the end of the data set. (It may be noted that 
in relation to nonstationary effects, the artefacts arising 
from the sample breaking leading to limited length of 
the time series is no different to the limited extent of 
time series due to other constraints as long as the data 
points closer to the break point are ignored.) The re- 
sulting stress-strain series for both the samples shown 
in Fig. and b are used for further analysis. 

We first investigate the dynamics of the less ductile 
Zr 65 Cui 5 Ni w Al w BMG. The total length of the time 
series used is about 29,000 points in units of St = 
0.02s. Before proceeding further, we first summarize 
the results of tests for stationarity (of the corrected time 
series), which involve calculating the mean, variance 
and autocorrelation function for several nonoverlapping 
subintervals of relatively long intervals. As the origi- 
nal time series is of finite length, the calculated quan- 
tities should match within some reasonable errors. We 
find these statistical quantities agree within a small error 
(for subintervals of length of about 5000 - 7000 points). 
Using only two intervals, each containing about 14500 
points, we have also calculated the autocorrelation func- 



tion for the first and second half of the time series of the 
data set. The calculated autocorrelation functions for 
the two halves were found to be nearly the same giving 
the same correlation time of 30 units as for the full data 
set (and therefore not shown). 

We have calculated the correlation integral, which is 
known to be more sensitive to nonstationarity Il36ll for 
the two halves of the data set and the full set. In gen- 
eral, using a much smaller delay time compared to the 
correlation time has been known to enhance the scaling 
regime. Thus, the delay time of 12 sampling units of 
time is used. A window size of w = 2(d - 1)t is used 
to exclude temporally correlated points [57]. A log-log 
plot of the correlation integrals for the first and second 
half is shown in Fig. |5]for t = 12 and for embedding 
dimensions d = 3 — 6. As can be seen, the slopes for 
successive embedding dimensions converge to ~ 2.2 for 
both halves of the data set. This turns out, as we shall 
see, to be nearly the same as the converged value of the 
slope for the full data set. Thus the "corrected" data set 
is taken to be stationary. 

To estimate the correlation dimension, we shall use 
the full data set. A plot of In C(r) verses In r for the 
full data set of 29, 000 points is shown in Fig. [6^ for 
embedding dimensions d = 3 — 6. It is clear that the 
slopes of In C(r) verses In r converge to about 2.2 as the 
embedding dimension reaches d — 6. Comparing these 
plots with those in Fig. [5] (for the two halves of the 
data set), we see that the local slope fluctuations in Fig. 
[6^ is much less than in Fig. [5] This feature is to be 
expected for a longer time series. The converged value 
of the slopes of In C(r) is taken to be the correlation 
dimension v ~ 2.2 of the signal. 

The Lyapunov spectrum is calculated (for the full 
data set) keeping d — 4 with a time delay of t — 2 
and a shell size 0.3%. The spectrum is shown in Fig. 
|6p. Note the existence of a stable positive and a zero 
Lyapunov exponent. The Kaplan- Yorke dimension cal- 
culated from the Lyapunov spectrum turns out to be 
Dk y ~ 2.84, which is consistent with the correlation 
dimension. The existence of a finite correlation dimen- 
sion and a stable positive (and zero) Lyapunov exponent 
strongly suggests that the stress signal is chaotic. (We 
have also calculated the Lyapunov spectrum for the two 
halves of the data set. The values of the positive and 
zero exponents are close to those of the full data set.) 

Now consider the dynamic analysis of the time series 
for the ductile CM47.5Zr47.5AZ5 BMG sample S d- The to- 
tal length of the time series used is about 18, 000 points 
in units of 8t = 0.05 s, which is quite short. While the 
plastic strain is much larger than the less ductile sam- 
ple S b, the coarse sampling rate imposes limitation on 
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Figure 5: (a and b) Log-log plot of the correlation integral as a func- 
tion of r for the the first and the second half of the stress-strain series 
for the less ductile Zr^CuisMioA/io. The delay time is 12. 
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Figure 6: (a) Log-log plot of the correlation integral as a function 
of r for the stress-strain series for the less ductile Zr^CuisNho^ho 
sample for d = 3 - 6 with delay time 12. (b) Plot of the Lyapunov 
spectrum for the same time series with d = 4 and delay time 2. 



the analysis. The auto-correlation time for the stress- 
time series is about 18 units of 5t. In this case, sys- 
tematic changes in the mean stress level have been re- 
moved by subtracting the running average and the am- 
plitudes are rescaled using a local average. The result- 
ing time series is expected to be stationary. The mean, 
variance of subintervals of the data set does not show 
much dispersion. The autocorrelation functions for the 
first and the second half of the time series remain nearly 
the same, giving the same correlation time as that for 
the full data set. This suggests that the time series is 
stationary. However, the stronger test for stationarity 
that should be reflected in the calculated correlation in- 
tegrals from the first and second half of the data set as 
with the full data set do not offer any additional sup- 
port as the slopes of In C(r) verses In r do not converge 
even though the slopes for each embedding dimension 
are nearly the same (not shown). In view of this, we 
present the results for the full data set that shows no ev- 
idence of chaos. We have used a range of typical time 
delay values along with suitable window sizes to calcu- 
late the correlation integral. However, we have not been 
able to find any convergence in the slopes of the corre- 



lation integral. A plot of In C(r) verses In r is shown in 
Fig. |7]for embedding dimensions d = 3 — 7. It is clear 
that there is no evidence of convergence of the slopes 
even for d — 7. A calculation of Lyapunov spectrum 
also did not show the existence of positive exponent or 
a zero exponent. On the basis of this we conclude that 
there is no evidence of chaotic dynamics. 

4.2. Self-organized critical dynamics 

A statistical analysis of the stress-strain curves for 
both the samples has also been been carried out. First, 
consider the less ductile sample S b- By visual inspec- 
tion, it is clear that most stress drops are large with much 
fewer small stress drops. This is what is reflected in 
the distribution of stress drop magnitude shown in Fig. 
[8^. The figure shows a nearly symmetric distribution 
around an average value of 3QMPa, but there is a peak 
for small values (about 2MPa). The corresponding dis- 
tribution of the durations also exhibits a similar form. 
No scaling regime could be detected in these distribu- 
tions. 
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Figure 7: (a) Log-log plot of the correlation integral as a function of r 
for the stress-strain series for the ductile Cw47.5Zf47.5Ai5 BMG sample 
for d = 3 - 7 with delay time 12. 



In contrast to the less ductile sample S b, a visual in- 
spection of the stress-time series for the ductile sample 
S d suggests that there are a large number of small stress 
drops with fewer large drops. This feature is similar to 
the hig h strain-rate type A serrations of the PLC effect 
ifsL 1 1 Oil wherein the analysis showed a power-law distri- 
bution of stress drops. Following this, we have calcu- 
lated the distribution of the stress drops as a function of 
the magnitude of the stress drops. As can be seen from 
Fig. |8j), the distribution of stress drops shows a mono- 
tonic ally decreasing trend, which when plotted on a log- 
log scale shows that the distribution follows a power law 
of the form D(Acr) ~ Acr~ a . This is shown in Fig. [9^. 
The scaling regime is nearly two orders in Act with an 
exponent value a ~ 1.5. 

Since the stress drop magnitudes are related to their 
durations, a similar analysis for the distribution of the 
durations of the stress drops At can be investigated for 
the ductile case. The distribution shows a tendency to 
follow a power-law scaling of the form D(At) ~ Ar 13 , as 
reflected in Fig. |9j). The exponent/? ~ 3.26. However, 
the extent of scaling regime is much less than that for the 
stress drop magnitudes. Such a limited scaling regime 
in the event durations is well documented. Indeed, even 
in model systems, such as the sand-pile model ||4^|6J], 
the extent of scaling regime in durations is much less 
than that for the magnitudes. 0. For model systems, it is 



3 For instance, in the large-scale simulations of sand-pile model 
|64|| . while the scaling regime for the magnitudes is nearly five order, 
that of durations is at best three orders. Further, in most experimental 
work that we are aware of where scaling for magnitudes and durations 
are reported, the scaling for duration is seldom more than one and half 
orders [50]. Quite often only the distribution of event sizes is taken to 
identify the process as SOC 165ll . 
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Figure 8: (a and b) Distribution of stress drop magnitudes for the 
less ductile Zr^Cii[sNiioAI[o and more ductile Cn47.5Zr47.5AZ5 BMG 
samples, respectively. 



easy to trace this features to the fact that the life time of 
a cluster is much shorter than its size 049II . In this case, 
since /3 (the exponent for durations) is high, sampling 
rate should be high for the same scaling extent to be 
realized compared to a situation where /3 ~ 1 as in the 
case of the PLC effect. In addition, in our case, part 
of contribution to limited scaling regime arises due to 
the coarse sampling rate (2QHz) and due to the limited 
precision in the measurement of the stress. We have 
also plotted the mean durations of the stress drops as a 
function of the size of the stress drops shown in Fig. \10\ 
The exponent value is x ~ 0.32. The scaling regime 
is nearly two orders as for the stress drop magnitudes. 
It can be easily checked that the scaling relation a = 
x(f3 - 1) + 1 is satisfied reasonably well. 

An independent check of the scaling relation can be 
carried out by computing the exponent corresponding to 
the power spectrum of the stress-strain curve. Kertesz 
and Kiss B66I1 derived an alternate relation between the 
exponents a and x that relates to the low-frequency be- 
havior of the power spectrum. By assuming that the 
total energy dissipated stems from independent elemen- 
tary events whose energy density spectrum is quasi- 
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Lorentzian, these authors have shown that if the scaling 
exponent satisfy the inequality 



-0.5 



2x + a > 3, 



(9) 



then the lower-frequency component of the power spec- 
trum S (of) should behave as: 



S((L>) ~ (jJ 



<3-a)/x 



(10) 



Thus, if 2x + a < 3, the S(u>) will always scale with 
l/co 2 . 

The power spectrum of the time series for the duc- 
tile alloy for the two halves and the full data set have 
been calculated. The low-frequency region for these 
three data sets can be well fitted to a power law with 
nearly the same exponent value ~ 1.9. Since the plots 
are practically the same, in Fig. QT|we show the power 
spectrum for the full data set. This is clearly consistent 
with the expected slope of l/a> 2 as 2x + a < 3. 

Despite the limited scaling regime in the time dura- 
tions, the cross check increases the confidence level in 
asserting the SOC nature of the stress-time series of the 
ductile sample. 
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Figure 9: (a and b) Distribution of stress drop magnitudes and dura- 
tions of the ductile C«47 iZr$-] 5A/5 sample, respectively. 
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Figure 1 1 : Power spectrum for the first half and the full data of the un- 
cleaned stress-strain series respectively of the ductile CM47 sZr^ 5A/5 
BMG sample. 



Finally, some comments may be useful on issues 
of nonstationarity. We note that the stress-time series 
of the ductile sample shows hardening, while that of 
the less ductile sample is nearly constant. However, 
the ductile data, which requires subtracting the run- 
ning average from the raw data set and rescaling, is not 
chaotic, but displays power law distribution. However, 
the power-law nature of the distributions remains the 
same for the detrended data set or the raw data set as 
the definition of the event is not influenced by drift of 
the mean stress level. Thus, the results of the analysis 
are not affected by subtracting the mean drift in stress- 
time series. In contrast, the less ductile sample data is 
nearly flat, but is chaotic, requiring stationary time se- 
ries. In this sense, removing the mean stress level from 
the time series has very little effect on our results. 
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5. Summary and discussion 

In summary, the analysis of the stress-strain curves 
from the intermittent flow of the two BMG samples 
shows that the deformation dynamics is very complex. 
Two complex but distinct dynamical regimes, namely 
chaos and SOC, are discovered. The underlying dynam- 
ics responsible for irregular serrations in the less duc- 
tile Zr^Cu\sNi\oAlu) BMG is chaotic; chaos requires 
a small number of collective degrees of freedom for a 
proper description of the dynamics. Indeed, since the 
associated correlation dimension is v ~ 2.2, the mini- 
mum number of degrees of freedom, i.e., the minimum 
number of variables required to describe the dynamics 
of serrations is given by the integer larger than v + 1, 
which in this case is four lf58ll . In contrast, for the duc- 
tile CM47.5Zr47.5AZ5 BMG, power law distributions are 
seen for stress drop magnitudes and durations. While 
two orders of magnitude for the extent of scaling for the 
stress drop magnitudes (and the mean durations of stress 
drops) is reasonable for any experimental data, the ex- 
tent of scaling for durations is only one order. Consider- 
ing the coarse sampling rate as well as the limited con- 
trol in measuring the stress, coupled with the fact that 
the exponents also satisfy the scaling relation and the in- 
dependent check using power spectrum, it appears that 
the ductile CM47.5Zr47.5AZ5 BMG follows SOC type of 
dynamics. Due to the scale-invariant nature of the pro- 
cess, SOC dynamics requires large number of degrees 
of freedom (actually infinite in theory), which in turn 
suggests that a hierarchy of length scales is responsi- 
ble for the stress-strain curves. Noting that stress-strain 
curves are a result of spatial average over all the shear 
band activity in the entire sample, the fact that two dif- 
ferent kinds of dynamics could be detected comes as a 
surprise. Nevertheless, the nature of the serrations ap- 
pears to contain this hidden information. The next step 
is to interpret these results in the light of available phys- 
ical information and models to see if further insight is 
possible. In the following we attempt to build a phys- 
ical picture that is consistent with the present state of 
knowledge. 

To this end, we note that there is considerable similar- 
ity between the serrated flow of the BMGs and the PLC 
effect even though the underlying physical processes are 
different. Thus, much can be inferred from the similari- 
ties and the differences between the serrated flow in the 
BMGs and the PLC effect H \ul Q . Recall that 
the PLC serrations and the associated bands arise due to 
collective pinning and unpinning of dislocations from 
solute atmosphere, corresponding to the stick and the 
slip phases, respectively. Thus, longer the waiting time, 



larger is the stress required to unpin them. For low strain 
rates, there is adequate time for aging to be completed, 
which therefore results in large yield drops. Increase in 
applied strain rate restricts the aging process, leading to 
a decrease in the serration amplitudes 

In the case of the PLC effect, at low strains, large am- 
plitude type C serrations are seen [J3 , 39 1 . Time-series 
analyses of the PLC type C serrations (in different sam- 
ple types) have revealed that the type C serrations are 
chaotic JESS]. On the other hand, the type C ser- 
rations are identified with randomly nucleating disloca- 
tion bands where each stress drop can be identified with 
a single band, i.e., the entire plastic strain rate is lo- 
calized to a small spatial extent lIM l39ll . This picture is 
substantiated by the Ananthakrishna model for the PLC 
effect The similarity with the large-amplitude 

chaotic serrations of the Zr^Cu^NiioAlio BMG sam- 
ple is clear. The fact that the PLC type C serrations are 
due to localization of dislocations into bands suggests 
that each stress drop in the case of the Zr^Cu^NimAlio 
BMG sample may be a result of localization of the shear 
to a small spatial region. This inference appears to be 
supported by recent observations on samples that ex- 
hibit significant ductility despite a single shear opera- 
tion. In such ductile BMGs, each stick-slip event has 
been shown to correspond to localized striations (of reg- 
ular spacing) II37[ l38l 15211 . These studies attribute the 
significant ductility to repetitive sliding of a major shear 
band. In contrast, a vein-like structure is usually found 
in low ductility situations where catastrophic failure is 
triggered by a single shear band. 

Indeed, Cheng et ai, have proposed a stick-slip 
model that is applicable to the instability process when a 
major shear band is observed [37]. These authors show 
that in this case the temperature rise is small, with the 
deformed sample remaining cold. While it is not easy 
to measure the in situ temperature inside the sample, an 
estimate based on their model turns out to be small for 
our case. Noting that the composition of our sample is 
almost the same as used in their experiment, it appears 
that their model adequately captures the results of the 
dynamics of serrations in the less ductile BMG. 

In general, stick-slip dynamics is considered to be 
regular, meaning that they can be approximated by peri- 
odic recurrence of stress drops (see Ananthakrishna and 
De in Ref. i67ll .) The present analysis, however, shows 
that the serrations are not regular but chaotic with mini- 
mum of four variables required for dynamic description 
of serrated flow ll58ll . Identifying these four degrees of 
freedom with physical variables should give some in- 
sight into the underlying dynamics. This will be done 
based on known physical mechanisms and model stud- 
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ies used in describing the deformation process. 

Several theoretical studies has demonstrated the im- 
portance of free volume and local temperature in under- 
standing the deformation of BMGs lEl I24H2I [26l [32 . 
3711 . It is well known that shear localization is affected 
by rise in the local temperature and volume dilation due 
to free volume generation. The local plastic strain rate 
depends on the available free volume, temperature and 
the local stress within the shear band. The latter is con- 
trolled by applied strain rate through the machine equa- 
tion (stress rate is proportional to the difference between 
the applied and plastic strain rate). Thus, the three obvi- 
ous variables that can be identified are the local plastic 
strain rate, free volume and temperature. However, both 
temperature and free volume in the shear band are influ- 
enced by the applied strain rate, which also affects the 
instantaneous flow stress in the shear band. Moreover, 
since the applied strain rate limits the extent of contri- 
bution from these three variables, which in turn depends 
on stress, the instantaneous stress in the shear band can 
be identified as the fourth variable. Thus, we rationalize 
that a minimal model should include local plastic strain 
rate, free volume, temperature and stress. Indeed, these 
are the four variables used in the model proposed by 
Cheng et al. 13711 . However, if the intended goal is to 
model the fluctuating stress-strain curves as the spatial 
average over the shear band activity in the entire sam- 
ple, this objective appears to be beyond the current state 
of knowledge. 

In the case of the ductile Cw47.5Zr47.5AZ5 BMG, the 
serrations are quite irregular with large number of small 
stress drops. This feature is similar to the PLC type 
A serrations seen at high strain rates for the propagat- 
ing bands |H H B EM 0. The small magnitude 
stress drops are attributed to incomplete plastic relax- 
ation. In the PLC type A serrations, power-law distribu- 
tions of stress drop magnitudes and durations have been 
reported it H OH] . Since the serrations from the ductile 
BMG also have a large number of small stress drops as 
in the PLC type A serrations, the possibility of a power- 
law SOC state is high. While the extent of the scaling 
is somewhat limited, an SOC type of dynamics is sug- 
gested. We note here that even in the case of PLC type 
A serrations, the sampling rate is not significantly better 
but the control on stress was improved. Given these lim- 
itations, the interpretation of the physical origin of the 
SOC-type dynamics could be somewhat speculative. 

To understand the origin of the power-law distribu- 
tion of stress drop magnitudes and durations in the duc- 
tile Cw47.5Zr47.5AZ5 BMG sample, consider Fig. [3j). In 
this case, multiple shear bands could be observed on the 
surface of the deformed sample. Further, the magnified 



view of a severely deformed region shown in Fig. |3j; 
clearly suggests a complex interaction among the large 
number of shear bands. This also suggests that the pos- 
sibility of a hierarchy of length scales over which the 
shear can propagate. Thus, a large number of metastable 
states can be created due to the interaction between the 
shear bands. This includes short-range interaction from 
the intersection of the bands and consequent arrest, and 
long-range interaction through the strain fields. Then, 
each stress drop corresponds to the system surmounting 
the barrier (corresponding to the arrested shear bands) 
jumping to a neighboring metastable state only to be 
caught in another metastable state. This picture is con- 
sistent with SOC-type dynamics. 

Liu et al. 15 311 have shown that the large plasticity 
is caused by the simultaneous nucleation of numerous 
shear bands (in different directions) and their interac- 
tions. The microstmcture of the ductile CM47.5Zr47.5AZ5 
BMG has been shown to consist of hard regions sur- 
rounded by soft regions. The wavelength of the wavy 
shear bands is consistent with the dimensions of the 
hard regions. They ascribe the plasticity of the ductile 
BMGs to concurrent nucleation and hindered propaga- 
tion of the shear bands throughout the sample. From 
this point of view, SOC-type of dynamics is likely to 
be seen in many other ductile metallic glasses [fl6i I2H 
El |H |H S ED, III Indeed, Wang et al Hshow 
that there is a scaling regime for the distribution of elas- 
tic energy density for several samples that are relatively 
more ductile. 

However, one basic difference between the stress- 
strain curves of the ductile BMG and the PLC type A 
serrations is the sharpness of the stress drops in the for- 
mer case as can be seen from Fig. |2j). Indeed, the 
reloading part is nearly linear except for a possible small 
curvature seen just before the stress drop. This must be 
contrasted with the stress-strain curves for the type A 
serrations of the PLC effect, where there is a substantial 
deviation from linearity, suggesting a plastic contribu- 
tion on the loading part. More importantly, there are 
small stress fluctuations both on the ascending part as 
well as the descending part of the stress-time curves of 
the PLC type A serrations. (See Fig. le of Ref. fl.) 
The latter feature is the signature of partial relaxation 
of plastic strain bursts. Thus, the sharp stress drop for 
the ductile BMG also implies that all shear band activity 
is close to the fully relaxed state. 



4 The basic idea used in designing ductile BMGs is to nucleate 
large number of shear bands that can interact with each other, thereby 
preventing a catastrophic failure. This can be accomplished, for exam- 
ple, by providing a heterogeneous microstructure such as nanoscale 
crystallinity |69, 73. IzH |72[| - 
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Finally, we note that while the less ductile sample is 
chaotic it does not follow SOC dynamics (i.e. the stress 
drop magnitudes and durations do not follow power- 
law distribution). In contrast, the ductile sample that 
follows SOC dynamics does not show any evidence of 
chaos. This is consistent with the well-established fact 
that these two dynamic states are considered to be very 
distinct: chaos is characterized by a few degrees of free- 
dom, while SOC is characterized by a large number of 
degrees (theoretically infinite) of freedom. Due to this 
fundamental distinction, these two types of dynamics 
are considered mutually exclusive and are not realized 
in the same system. Indeed, in the PLC effect also, the 
chaotic type C and B serrations are observed at low and 
medium strain rates, while type A serrations that fol- 
lows SOC dynamics are observed at high strain rates 
USE. 

In summary, we have analyzed the serrated flow be- 
havior for the two typical BMGs by time-series analysis 
and statistical methods. For the less ductile alloy, we 
have demonstrated the chaotic nature of the stress serra- 
tions by showing the existence of a finite correlation di- 
mension and a positive Lyapunov exponent. In contrast, 
for the ductile alloy, the distributions of stress drop mag- 
nitudes and their time durations appear to follow power- 
law scaling form reminiscent of self-organized critical 
state. The three exponents also satisfy the scaling rela- 
tion that should be satisfied by power-law SOC dynam- 
ics 1 49[ 51]. To the best of our knowledge, this is the first 
time that the two dynamical regimes are observed in the 
deformation of the two different BMGs. Considering the 
fact that the extent of plasticity is only a factor of two 
between the two samples, and the fact that the method 
of analysis uses only the stress-strain curves, this dif- 
ference in the nature of dynamics is surprising. Indeed, 
that the method is able to detect this difference is also 
surprising. However, we have suggested that the ba- 
sic difference in the dynamics can be related to the dif- 
ferences in the physical mechanisms contributing to the 
different ductility levels in the two BMGs. The chaotic 
dynamics of the less ductile BMG appears to be consis- 
tent with the accepted picture of a single shear band slid- 
ing in repetitive manner, leading to considerable ductil- 
ity. Each stick-slip corresponds to localized striations of 
regular spacing in samples 11371 B8l 15211 . Further, for this 
case, the minimum number of variables required for the 
description of the chaotic serrations is four. On the ba- 
sis of known physical mechanisms, we have identified 
these variables to be the free volume, temperature, local 
plastic strain rate and stress; this is consistent with the 
model of Cheng et al 13711 . In contrast, for the ductile 
BMG, the high ductility appears to be related to concur- 



rent nucleation of large number of shear bands through- 
out the sample, which in turn can give raise to a hier- 
archy of length scales over which shear bands propa- 
gate. Our analysis also suggests that the ductile nature 
emerges from the production of several shear bands and 
their interactions rather than partial plastic relaxation of 
the shear activity in the shear bands as in the PLC effect. 
These rationalizations may have implication in under- 
standing the mechanism controlling plastic deformation 
of the ductile BMGs and should prove useful in model- 
ing efforts. 
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